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Grosman  B,  Shaik  OS,  Helwig  BG,  Leon  LR,  Doyle  FJ  3rd.  A 

physiological  systems  approach  to  modeling  and  resetting  of  mouse  thermo¬ 
regulation  under  heat  stress.  J  Appl  Physiol  111:  938-945,  2011.  First 
published  June  23,  2011;  doi:  10. 1 152/japplphysiol.005 19.2010. — Heat 
stroke  (HS)  is  a  serious  civilian  and  military  health  issue.  Due  to  the  limited 
amount  of  experimental  data  available  in  humans,  this  study  was  conducted 
on  a  mouse  mathematical  model  fitted  on  experimental  data  collected  from 
mice  under  HS  conditions,  with  the  assumption  there  is  good  agreement 
among  mammals.  Core  temperature  (Tc)  recovery  responses  in  a  mouse 
model  consist  of  hypothermia  and  delayed  fever  during  24  h  of  recovery  that 
represent  potential  biomarkers  of  HS  severity.  The  objective  of  this  study  was 
to  develop  a  simulation  model  of  mouse  Tc  responses  and  identify  optimal 
treatment  windows  for  HS  recovery  using  a  three-dimensional  predictive  heat 
transfer  simulation  model.  Several  bioenergetic  simulation  variables,  includ¬ 
ing  nonlinear  metabolic  heat  production  (W/m3),  temperature-dependent 
convective  heat  transfer  through  blood  mass  perfusion  (W/m3),  and  activity- 
related  changes  in  circadian  Tc  were  used  for  model  simulation.  The  simu¬ 
lation  results  predicted  the  experimental  data  with  few  disparities.  Using  this 
simulation  model,  we  tested  a  series  of  ambient  temperature  treatment 
strategies  to  minimize  hypothermia  and  delayed  fever  to  accelerate  HS 
recovery.  Using  a  genetic  algorithm,  we  identified  eight  time  segments 
(ambient  temperature  =  27,  30,  31,  29,  28,  28,  27,  26°C)  of  110  min  total 
duration  that  optimized  HS  recovery  in  our  model  simulation. 

numerical  modeling;  heat  stroke;  hyperthermia;  hypothermia;  genetic 
algorithm 


mammalian  species  have  developed  numerous  mechanisms  to 
cope  with  a  wide  variety  of  environmental  stressors,  including 
heat  stress.  Heat  stress,  which  can  progress  to  fatal  heat  stroke 
(HS)  is  a  serious  public  and  military  health  issue  (16).  HS 
results  from  prolonged  exposure  to  a  hot  environment,  causing 
core  temperature  (Tc)  to  exceed  the  normal  temperature,  lead¬ 
ing  to  severe  hyperthermia  that  can  result  in  physical  and 
cognitive  function  decrements.  As  such,  HS  is  clinically  char¬ 
acterized  by  hyperthermia  (Tc  typically,  but  not  always  greater 
than  40°C)  collapse  and  profound  encephalopathy  that  presents 
as  delirium,  agitation,  stupor,  seizures,  or  coma  (6).  Rapid 
cooling  therapy  and  fluid  resuscitation  represent  the  current 
treatment  strategies,  but  there  are  currently  no  effective  phar¬ 
macological  treatments  for  HS.  Despite  clinical  cooling  thera¬ 
pies  that  rapidly  lower  body  temperature  to  minimize  tissue 
injury,  HS  is  often  fatal  and  results  in  permanent  neurological 
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damage  in  30%  of  survivors  (3).  Identification  of  novel  thera¬ 
pies  and  optimal  treatment  windows  to  accelerate  recovery  are 
needed  to  minimize  HS  morbidity  and  mortality. 

In  the  present  study,  we  used  a  mathematical  model  to 
simulate  mouse  Tc  responses  during  heat  stress  and  recovery 
and  identify  optimal  time  windows  for  ambient  temperature 
(Ta)  treatment  strategies  to  optimize  resetting  of  the  thermo¬ 
regulatory  profiles  back  to  normal  Tc.  In  the  past  three  decades, 
several  mathematical  models  of  the  human  thermoregulation 
system  have  been  developed  to  enhance  understanding  of 
regulatory  processes  (5,  9,  4a,  18).  Most  of  these  models  were 
based  on  cylindrical  approximations  of  the  body  segments  with 
shell  layers  representing  the  body  tissue  layers.  Due  to  the 
limited  amount  of  experimental  data  available  in  humans  and 
the  inability  to  study  Tc  recovery  responses  in  HS  patients,  it 
has  not  been  possible  to  utilize  these  models  to  test  the  efficacy 
of  treatment  strategies  during  recovery.  In  this  study,  mouse  Tc 
responses  during  heat  exposure  and  40  h  of  recovery  were 
modeled  to  expand  current  prediction  capabilities  through 
several  circadian  cycles  and  test  the  efficacy  of  several  Ta 
treatments  on  HS  recovery. 

Our  laboratory  recently  developed  an  in  vivo  mouse  HS 
model  using  conscious  free-moving  animals  with  implanted 
radiotelemetry  devices  to  examine  Tc  responses  through  48  h 
of  HS  recovery  (14).  This  was  the  first  in  vivo  HS  model  to 
remotely  measure  Tc  through  several  circadian  cycles  without 
compromising  the  subject’s  natural  behavioral  and  autonomic 
thermo-effector  responses  to  prolonged  heat  exposure  (14). 
Following  heat  exposure  to  a  maximum  Tc  (Tc.max),  ranging 
from  42.4  to  42.7°C,  mice  exhibited  a  biphasic  thermoregula¬ 
tory  recovery  response  consisting  of  initial  hypothermia  fol¬ 
lowed  by  fever  at  24  h  (12,  13).  Survival  rates,  depths  of 
hypothermia,  and  recovery  times  of  mice  subjected  to  heat 
stress  were  mainly  affected  by  Tc,max  attained  during  heat 
exposure,  as  well  as  recovery  Ta. 

MATERIALS  AND  METHODS 

Mouse  heat  transfer  mechanisms  were  modeled  with  consideration 
of  conduction,  convection,  and  radiation  through  the  skin  to  the 
environment  as  a  passive  system.  In  addition,  metabolic  heat  produc¬ 
tion,  blood  flow,  saliva  spreading,  and  shivering  (referred  to  as  the 
active  controlled  system)  were  addressed,  in  part,  in  our  model.  The 
passive  system  was  modeled  as  multiple  layers  of  tissue  with  distinct 
anatomical  differences,  and  the  active  system  was  designed  to  be  a 
function  of  the  difference  among  multiple  tissue  temperatures  and  the 
desired  set-point  temperature  (T.ret).  In  homeothermic  animals,  main- 
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tenance  of  a  relatively  constant  Tc  is  achieved  by  an  increase  in 
metabolic  heat  production  and  reduction  in  heat  loss  when  Tc  <  Tset, 
and  a  reduction  in  metabolic  heat  production  and  an  increase  in  heat 
loss  when  Tc  >  T.set  (21).  The  difference  between  the  Tc  of  tissues  and 
Tret  serves  as  the  driving  force  for  heat  transfer  mechanisms. 

Mouse  geometry.  Mouse  geometry  was  idealized  as  an  ellipsoid  for 
the  torso  (length,  0.08  m,  diameter,  0.003  m)  and  a  cone  for  the  tail 
(length  0.075  m),  with  multiple  tissue  layers  (skin,  fat,  muscle,  core, 
and  bone).  The  thickness  of  tissue  layers  (skin,  fat,  muscle,  and  core 
organs)  was  based  on  approximate  percentage  of  the  total  body 
weight.  Metabolically  active  components,  including  the  brain,  heart, 
liver,  and  kidneys,  were  considered  in  the  core  layer  of  the  mouse 
torso.  Trakic  et  al.  (22)  yielded  a  25.7-g  mouse  voxel  phantom  with  a 
resolution  of  0.13  X  0.13  X  0.143  mm3  from  similar  0.39  X  0.39  X 
0.43  mm3  voxel  phantom  of  a  male  Sprague-Dawley  rat  with  a  mass 
of  ~370  g.  We  used  this  approximation  with  volume  of  3.65,  3.97, 
4.41,  and  4.50  mm3  of  the  skin,  fat,  muscle,  and  stomach  layers, 
respectively.  Volumes  of  brain,  heart,  liver,  and  kidneys  were  calcu¬ 
lated  from  the  preceding  voxel  data,  and  equivalent  volumes  of 
spheres  were  considered  in  the  model.  The  software  package  Comsol 
Multiphysics  was  used  to  define  the  geometry  and  mesh  the  model. 
Input  variables  included  initial  conditions  of  the  mouse  Tc,  global 
expressions  (Comsol  types  of  expressions  that  are  used  by  different 
domains  of  the  model),  and  boundary  conditions,  along  with  material 
properties.  Comsol  uses  an  adaptive  finite-element-based  method  to 
solve  the  resulting  partial  and  ordinary  differential  equations.  To 
minimize  the  computational  costs,  a  symmetric  problem  with  one-half 
of  the  actual  geometry  was  implemented. 

The  model  included  the  major  heat  transfer  mechanisms  of  con¬ 
duction,  radiation,  blood  perfusion,  metabolic  heat  production,  evap¬ 
oration,  and  convection.  The  temperature  distribution  inside  the  tissue 
was  computed  by  the  Pennes  bio-heat  transfer  equation  (BHTE)  (17). 
The  differential  equation  that  describes  the  heat  dissipation  in  a 
homogenous,  infinite  tissue  volume  was  given  by  the  following 
Pennes  BHTE  (7): 

3T(£,f) 

P(S)Cp(D — — —  =  K(QV2 T(£,  t) 

+  B(£,  t)  +  g(£.  t)  +  QSL  0 

where  £  =  (jc,  y,  z),  p  (kg/m3)  is  the  specific  mass  density,  Cp 
(J-kg-1-°c-1)  is  the  specific  heat,  T  (°C)  is  the  temperature  at  time  r, 
K  (W-m_1-°C_1)  is  the  thermal  conductivity,  B  (W/m3)  is  the  con¬ 
vective  heat  exchange  between  blood  and  tissue,  Q  (W/m3)  is  the 
metabolic  heat  generation,  and  Qr  (W/m3)  is  thermal  radiation  gov¬ 
erned  by  the  Stefan-Boltzmann  law.  Pennes  postulated  that  the  total 
heat  exchange  between  the  blood  and  the  surrounding  tissue  could  be 
modeled  as  a  nondirectional  heat  source  whose  magnitude  is  propor¬ 
tional  to  the  volumetric  blood  flow  and  the  difference  between  the 
local  tissue  temperature  and  the  blood  temperature  (17).  Blood  enters 
each  tissue  compartment  at  the  temperature  of  the  central  blood 
compartment  and  returns  to  the  central  compartment  at  the  current 
tissue  temperature.  The  central  circulatory  system  is  conceptual  and 
actually  simulates  the  blood  system  as  a  single  compartment,  in  which 
blood  from  various  organs  is  mixed,  then  redistributed.  The  rate  of 
heat  transfer  by  blood  flow  through  different  tissue  layers  is  modeled 
as  a  response  to  the  physiological  demands  and  thermoregulation  (22). 

Blood  perfusion  rates  to  different  tissues  is  temperature  dependent 
at  the  space  time  point,  which  is  controlled  by  combinations  of  global 
signals  (depending  on  temperature  receptors  throughout  the  whole 
body,  integrated  in  the  hypothalamus)  and  local  signals.  The  blood 
perfusion  equilibrium  rate  is  assumed  to  double  for  each  4.5°C  rise  in 
temperature  T(ij,r)  (2).  The  equilibrium  is  assumed  to  be  reached 
exponentially  by  the  following  time-dependent  equation  (23): 


—f  = - {ui,oexP[a(T  -  T0)]  -  w*}  (2) 

®  *  ^blood 

where  cov  and  uix,o  (m3-kg~'-s  ')  are  dynamic  and  basal  blood 
perfusion  [estimated  from  Trakic  et  al.  (22)]  for  a  specific  organ-tissue 
x,  respectively;  Tbiood  is  time  constant  (12  min);  and  a  =  ln(2/ 
4.5)°C-‘. 

Metabolism.  Resting  energy  expenditure  (REE)  or  basal  metabolic 
rates  are  expressed  as  a  function  of  mammalian  body  mass,  according 
to  Kleiber’s  law  (24).  An  organ-tissue  level  REE  model  is  formulated 
for  metabolically  active  organs,  such  as  the  liver,  heart,  brain,  kidneys, 
and  remaining  tissues.  Among  all  mammals,  the  brain,  liver,  kidneys, 
and  heart  account  for  73%  of  total  REE  (4).  The  equations  for 
metabolically  active  organs  and  remaining  tissues  are  based  on  Wang 
et  al.  (24) 

x  Tset 

2baS,,,o  =  «,-M-P-2  '  10  (3) 

where  Qbas^.o  (W/m3)  is  the  basal  metabolic  heat  production  of  the 
organ-tissue  x,  M  (kg)  is  the  body  mass  of  the  mammal,  ax  and  [3*  are 
specific  constants  of  the  organ  x,  and  T,  is  the  organ-tissue  tempera¬ 
ture  (22).  In  endotherms,  the  dependence  of  biochemical  reactions  on 
temperature  is  described  by  the  van’t  Hoff  Qio  effect  (24)  with  a 
sensitivity  coefficient  of  2. 

In  addition  to  blood  perfusion  rates  (Eq.  2),  thermoregulatory 
responses  consisted  of  metabolic  heat  production  that  was  based  on 
integrating  the  error  between  Tc  and  Tset,  given  by  the  following 
equations: 

(2ba v  v(  ( ■  I)  (2x,bas,0  ’  T/ , 


dMri 

dr 

kJJset  ~  Tc) 

Mr  l 

-Mr,low 

<Mr<  M, up 

Mr  =  - 

^r,up 

Mrl  > 

^r,up 

low 

Mn< 

Mr.iow 

where  Mr i  is  the  metabolic  regulation  factor,  which  is  driven  by  the 
deviation  between  Tc  and  Tset,  Mr  is  the  constraint  metabolic  regula¬ 
tion  factor,  and  km  (°C~Ts_1)  is  an  empirical  gain.  The  Mr  value  is 
bounded  between  Mr, iow  and  Mr.up,  the  lower  and  upper  constraints  to 
the  metabolic  regulation  factor,  respectively.  The  values  of  km  can 
differ  from  organ  to  organ  and  for  different  tissues;  however,  for  the 
present  case,  due  to  lack  of  biological  data,  it  is  assumed  to  be  equal 
for  all  tissues  and  organs  in  the  mouse  body.  The  exact  mechanism  of 
the  metabolic  regulation  is  unknown,  but  we  assumed  that,  after  a 
prolonged  heat  load,  the  body  is  trying  to  minimize  metabolic  activity. 
On  the  other  hand,  when  the  body  temperature  is  low  for  a  prolonged 
time,  Mt  can  get  the  value  of  Afr>up.  This  resembled  a  constraint 
integrator  in  process  control.  Notice  that  Mr  is  a  multiplier  on  the 
temperature-dependent  metabolic  function,  as  shown  in  Eq.  5. 

The  metabolic  heat  production  of  the  muscle  tissue  is  given  by  the 
sum  of  the  basal  value  <2bas,v  and  the  additional  heat  A<2m,  which  may 
be  produced  by  local  autonomic  responses,  such  as  normal  activity 
and/or  shivering  thermogenesis: 

2x(£>  t)  =  Gbas.muscle  +  A2ra 

activity  (5) 

*£bas,muscle  *  JQQ 

where  the  change  in  basal  metabolism,  A<2m,  is  the  difference  between 
the  actual  basal  rate  and  the  basal  rate  corresponding  to  neutral 
thermal  conditions,  and  activity  is  measured  as  number  of  gross  motor 
movements  per  minute.  Factors,  including  the  rate  of  movement,  can 
affect  activity  readings;  thus  reported  activity  should  be  considered  a 
relative  value  and  does  not  provide  insight  regarding  the  absolute 
distance  the  animal  moved.  Activity  was  simultaneously  collected  for 
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each  animal  throughout  the  data  collection  period  and  reported  by  the 
data  analysis  software  on  a  per-minute  basis  following  collection  of  Tc 
values.  We  calculated  the  correlation  coefficient  between  activity  and 
oxygen  consumption  on  570  experimental  data  points  (12)  and  found 
it  to  be  0.54  (P  value  of  2e-46),  meaning  that  there  is  a  significant 
correlation  between  metabolic  rate  and  mouse  activity. 

Circadian  variation  in  Tc  is  modeled  by  considering  the  activity 
levels  of  mice  and  the  drinking  of  water  at  25°C  (25).  An  additional 
heat  term,  <2m.a,  is  added  to  the  muscle  layer  and  other  metabolically 
active  organs  to  account  for  the  increase  in  Tc  due  to  activity  levels. 
Heat  is  lost  in  the  core  layer  of  the  mouse  body  due  to  the  drinking  of 
25°C  water.  This  amount  of  heat  (Qm.w,  J/m3)  is  removed  from  the 
core  layer  of  the  mouse  body  based  on  following  equation: 

Qm,  w  Pw  '  Cp,H2o(T  average  ^amb)  '  (0) 

*  stomach 

where  pw  (kg/m3)  is  the  density  of  the  water  consumed  by  the 
mouse,  Cp,h20  (J-kg_1-°c)  is  the  specific  heat  of  water,  Vwater  is 
the  amount  of  water  consumed  (m3),  and  VstomaCh  is  the  volume  of 
the  stomach  (m3). 

Heat  is  exchanged  from  the  mouse  body  by  convection  with  the 
ambient  air  ( Qc ;  W/m2),  by  water  evaporation  due  to  saliva  spreading 
on  the  skin  or  fur  surface  (Qe\  W/m2),  and  by  respiratory  losses  (Qrp\ 
W/m2).  The  rate  of  heat  exchange  over  the  body  of  the  mouse  varies, 
with  regard  to  the  furred  torso  and  the  unfurred  appendages.  Hence, 
heat  balances  from  furred  torso  and  unfurred  tail  are  modeled  as 
different  boundary  conditions.  In  general,  the  net  heat  flux  Qs  (W/m2) 
exchanged  from  the  surface  of  the  body  is  equivalent  to  the  sum  of 
these  individual  heat  exchanges  ( Qs  =  Qc  +  Qe  +  Qip )■  The  heat  loss 
due  to  forced  convection  is  governed  by  the  equation: 

2C(U)  =  -MT(U)-Ta]  (7) 

where  <2c(£,t)  (W/m2)  is  the  convective  heat  loss,  h  (W-m_2-°C_1)  is 
the  convective  heat  transfer  coefficient,  Ta  is  in  °C.  In  this  work,  h  is 
evaluated  by  the  empirical  formula  given  by  Wooden  and  Walsberg 
(26).  Different  values  of  h  are  used  for  the  furred  torso  and  the  tail  ( h 
is  larger  by  25%  on  the  tail).  When  Tc  <  T.set,  mice  piloerect  their  fur 
to  create  an  extra  layer  of  insulation,  thereby  conserving  body  heat.  In 
contrast  when  Tc  >  T.set,  as  was  present  in  the  present  study,  mice 
relax  the  piloerection  response  to  facilitate  heat  loss  to  the  surround¬ 
ings.  This  behavioral  response  was  incorporated  into  the  model  by 
using  h  as  a  linear  function  of  average  temperature  on  the  torso  ( Eq . 
8)  that  results  in  h  —  3  W-m~2-°C~1,  when  Tc  =  Tset- 

h  =  (Tc  -  10)/9(W  ■  m_2-°C_1)  (8) 

In  rodents,  the  behavioral  spreading  of  saliva  on  the  ventral  body 
surfaces  is  an  essential  mechanism  for  heat  loss  through  water  evap¬ 
oration.  <2e  represents  evaporative  heat  transfer  rate  per  unit  area  and 
consists  of  passive  diffusion  and  evaporation  of  saliva.  Passive  diffu¬ 
sion  is  a  small  quantity  of  residual  mass  transfer  that  occurs  between 
skin  and  surroundings,  which  is  negligible  in  the  present  case.  Mouse 
Tc  increases  when  sufficient  heat  cannot  be  transferred  to  the  envi¬ 
ronment  by  the  normal  heat  loss  mechanisms  of  convection,  conduc¬ 
tion,  diffusion,  and  radiation.  The  result  is  an  increased  Tc  above  T.set, 
triggering  the  active  thermoregulatory  control  system  of  saliva  spread¬ 
ing.  When  exposed  to  heat  stress  near  40°C,  rodents  exhibit  a  peak  in 
evaporative  heat  loss;  however,  when  exposure  continues,  the  evap¬ 
orative  heat  loss  falls  rapidly  in  hamsters,  whereas  in  rats  it  is 
maintained  for  several  more  hours  before  heat  exhaustion  occurs  (7). 
Due  to  lack  of  data  available  for  mice,  the  salivation  profile  is  set  by 
optimization  to  best  fit  the  experimental  data,  and  it  is  triggered  by  the 
active  system  once  the  heating  begins.  Evaporative  heat  loss  is 
represented  by  the  following  differential  equation: 


k-H,o  d ms 

Qe  =  — - ! 

Asl  d  t 


(9) 


where  Asl  (m2)  is  the  salivary  evaporation  skin  surface;  ms  is  amount 
of  saliva  spreading;  and  \h20  =  2,256  X  103  J/kg  is  the  heat  of 
evaporation  of  water.  In  the  present  work,  only  the  front  ventral  part 
of  the  mouse  body  and  front  legs  are  considered  as  Asi.  Evaporative 
heat  loss  becomes  greater  as  the  gradient  temperature  between  Tc  and 
Ta  increases,  albeit  it  reduced,  after  the  dehydration  point  is  reached, 
eventually  becoming  ineffective  as  the  remaining  saliva  on  the  body 
surface  evaporates.  The  time  corresponding  to  the  point  of  dehydra¬ 
tion  (fd;  s),  duration  (p,;  s),  and  the  maximum  amount  of  saliva 
spreading  (ms<max;  kg/s)  at  any  time  during  heat  stress  are  parameters 
that  differ  from  mouse  to  mouse.  The  exponential  expression  that 
represents  the  heat  loss  due  to  saliva  spreading  is  described  by  Eq.  10. 

\ 

Qe  ~7  maxC 
4  si 


”d\2 

I1  , 


(10) 


Respiratory  heat  losses  occur  when  air  is  inhaled  at  Ta  and  exhaled 
close  to  Tc  with  high  humidity.  In  Schmidt-Nielsen  et  al.  (20),  it  is 
discussed  that,  in  small  rodents  (Kengero  rat),  the  mechanism  of  heat 
loss  through  respiration  is  highly  efficient,  and,  therefore,  it  can  be 
approximated  to  be  zero.  Therefore,  in  the  present  model,  we  do  not 
define  the  pulmonary  tract  as  a  separate  organ,  and  the  respiratory  heat 
losses  are  neglected. 


RESULTS 

The  mathematical  model  in  the  preliminary  stage  needs  to  be 
calibrated  to  the  experimental  data.  The  following  parameters 
were  left  as  degrees  of  freedom  and  were  determined  with  the 
help  of  optimization  [genetic  algorithm  (GA)]:  the  time  corre¬ 
sponding  to  the  point  of  dehydration  (fd;  s),  the  salivation 
duration  (|x;  s),  the  maximum  amount  of  saliva  spreading 
(iws,max;  kg/s1),  Tvet  in  the  experimental  starting  time,  the 
lowest  value  of  metabolic  regulation  factor  (Mr, iow),  the  em¬ 
piric  gain  km  (Eq.  4),  and  Mr o  (the  initial  value  of  Eq.  4). 
Table  1  specifies  results  of  the  optimization  and  the  parametric 
range  in  which  the  optimization  search  was  done. 

Tset  was  selected  as  a  parameter  in  the  calibration  of  the 
model  because  the  mice  showed  decreasing  Tc  before  heat 
stress  [i.e.,  during  transition  into  the  lights-on  (inactive)  pe¬ 
riod],  which  was  explicable  as  a  change  in  Tve,.  The  parameters 
of  the  salivation  (fdj  |x,  and  mSJllax)  were  chosen  because  of  its 
large  effect  on  mouse  thermoregulation  and  the  lack  of  exact 
measurements  in  this  species.  The  salivation  profile  using  the 
parameters  that  are  presented  in  Table  1  predicts  an  overall 
salivation  of  2.9  g  over  the  4  h  of  heating.  The  parameters  of 
Eq.  4  were  calibrated  due  to  a  lack  of  experimental  data 


Table  1.  List  of  parameters  used  to  calibrate  the 
mathematical  model  to  the  experimental  data 


Parameter 

Range 

Value 

fd,  s 

2,400-6,000 

4,425 

|Jl,  s 

6,667-20,000 

1.0667e  +  004 

ms, max,  kg/s 

1.33 12e-007-3.3280e-007 

2.4545e-007 

Ts„,  °C 

32-38 

35 

^fr.low 

1-2 

1.625 

km,  °C~1*S~1 

le-5-4e-5 

2.875e-5 

Mt o 

2-3 

2.6875 

Shown  are  the  parameter  denomination,  range  of  optimization  search,  and 
final  optimization  value.  See  text  for  definition  of  parameters. 
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Fig.  1.  The  mathematical  model  response  vs. 
two  experimental  data.  The  figure  depicts  the 
good  agreement  of  the  mathematical  model 
(solid  line)  to  two  curves  of  empirical  data 
(dotted  and  dashed  lines).  To  calibrate  the 
mathematical  model  to  the  experimental  data, 
seven  parameters  were  adapted  using  optimi¬ 
zation.  The  typical  behavior  that  is  observed 
in  the  stage  of  the  heating  is  controlled  by 
regime  of  the  salivation  in  combination  with 
set-point  temperature  (Tvet)  that  is  fixed  for 
34°C.  Tc,  core  temperature. 


describing  the  function  that  regulates  metabolic  activity  in  the 
presented  mathematical  model;  however,  the  value  of  km  can 
be  approximated  by  the  speed  of  the  regulation  that  was 
observed  in  the  experiments. 

Figure  1  shows  the  successful  agreement  between  the  cali¬ 
brated  mathematical  model  and  the  experimental  data.  The 
variation  of  Tc  using  Tvel  =  34°C  from  the  beginning  of  the 
simulation  until  the  mouse  is  removed  from  the  heat  chamber 
(4  h),  and  Tvet  =  36.35°C  for  the  rest  of  the  simulation,  results 
in  a  very  good  agreement  with  the  experimental  data.  The 
constrained  metabolic  regulation  factor  (Mr i)  dynamic  behav¬ 
ior  is  depicted  in  Fig.  2.  We  predicted  that,  in  order  for  the 
mouse  to  achieve  homeostasis,  the  metabolic  heat  production 
that  is  predicted  by  Kleiber’s  law  (known  to  change  between 
different  species)  should  be  multiplied  by  approximately  three. 
Moreover,  under  heat  stress,  the  mouse  was  predicted  to  reduce 
its  metabolic  heat  production  by  ~50%.  We  predict  that  the 
main  reason  hypothermia  occurred  during  recovery  was  due  to 
a  reduction  in  metabolic  activity,  which  did  not  show  a  reversal 
for  several  hours. 

A  coarse  mesh  of  the  mouse  model  can  be  seen  in  Fig.  3.  As 
a  first  case,  a  passive  cooling  scenario  was  simulated  without 
considering  the  additional  heat  production  and  blood  perfusion 
terms  in  the  overall  BHTEs.  At  the  end  of  HS  experimentation, 


a  subset  of  mice  (N  =  4)  was  killed  and  placed  into  a  preheated 
environmental  chamber  until  a  Tc  of  42.7°C  was  attained. 
Following  this  heating  paradigm,  mice  passively  cooled  at  Ta 
of  25°C  for  60  min.  Simulation  results  of  the  computational 
model  are  compared  with  these  experimental  passive  cooling 
data  in  Fig.  4.  Tc  profiles  of  the  simulated  model  reached  Tc  = 
32°C  from  Tc  =  42.7°C  slowly  compared  with  experimental 
data  for  passive  cooling.  Sensitivity  analysis  of  the  model  with 
respect  to  parameter  h  was  performed  by  varying  ±20%  of  its 
nominal  value  of  h  =  3.5  W  m_2-°C-1.  By  increasing  and 
decreasing  the  value  of  h,  passive  cooling  times  were  de¬ 
creased  and  increased,  respectively;  however,  the  deviation  of 
Tc  from  the  nominal  cooling  curve  was  not  significant.  Sensi¬ 
tivity  analysis  of  the  model  was  also  carried  out  with  respect  to 
scaling  of  the  mouse  body.  Whole  mouse  geometry  was  scaled 
by  ±10%  from  its  nominal  size;  mice  with  10%  reduction  in 
scale  showed  similar  passive  cooling  profiles  to  those  from  the 
experiments.  The  nominal  mouse  geometry  was  updated  with  a 
10%  reduction  in  size,  and  the  new  dimensions  were  used  for 
further  calculations. 

The  spatial  temperature  distribution  across  a  mouse  body  for 
the  hyperthermia  region  at  t  =  200  min  is  shown  in  Fig.  5. 
When  exposed  to  Ta  =  39.5°C,  the  tail  temperature  rose 
quickly  due  to  increased  skin  blood  flow  to  a  heat-dissipating 


Fig.  2.  Variability  of  the  constrained  meta¬ 
bolic  regulation  factor  (Mr)  during  stage  of  the 
preheating,  the  heating,  and  the  acclimatiza¬ 
tion  periods.  During  the  heating  and  in  the 
early  stage  of  the  acclimatization,  Mr  reach 
the  lower  bound;  later  it  takes  Mr  around  3  h 
to  arrive  to  its  upper  bound,  while  recovering 
from  the  hypothermia. 
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Fig.  3.  A  three-dimensional  computational  mouse  geometry  with  torso  and  tail 
is  shown.  Dimensions  of  the  torso  and  tail  are  taken  from  average  experimental 
measurements  of  mice  subjected  to  heat  stress.  The  torso  consists  of  four 
layers:  skin,  fat,  muscle,  and  core.  The  tail  layer  consists  of  core  and  skin. 
Metabolically  active  organs,  heart,  liver,  brain,  and  kidney,  are  considered 
spheres  in  core  layer  of  the  torso.  Only  one-half  of  the  mouse  geometry  is  used 
in  the  computational  domain  due  to  symmetry  of  the  mouse  body  around  the 
center  and  to  minimize  computational  cost.  A  finite-element  mesh  with  44,267 
nodes  is  generated  for  the  computational  mouse  model  in  Comsol. 


Max:  38.304 


Min:  15.179 


organ  with  a  larger  surface  area-to-volume  ratio  and  higher  h 
value.  Heat  loss  due  to  spreading  of  saliva  on  the  front  ventral 
part  and  feet  of  the  mouse  body  were  observed  in  the  surface 
plots.  A  maximum  temperature  gradient  of  23°C  was  observed 
from  the  center  of  the  body  surface  to  front  legs  due  to 
evaporative  heat  loss  from  spreading  of  saliva.  On  average,  the 
front  part  of  the  mouse  body  was  exposed  to  lower  tempera¬ 
tures  compared  with  the  tail  and  center  part  of  the  body,  where 
liver,  kidneys,  and  intestines  reside.  Within  the  hyperthermia 
region,  temperatures  of  these  organs  were  more  than  the  actual 
mouse  Tc,  suggesting  these  internal  organs  were  more  suscep¬ 
tible  to  thermal  injury.  Also,  due  to  the  high  metabolic  activity 
in  the  liver,  the  average  liver  temperature  is  assumed  to  be 
more  than  the  mouse  Tc  observed  in  simulations. 

The  three-dimensional  thermoregulation  model  was  used  to 
optimize  the  recovery  time  of  a  heat-stressed  mouse,  by  vary¬ 
ing  Ta  in  the  hypothermia  region.  Ta  was  optimized  in  the 
control  region  consisting  of  eight  110- min  segments,  totaling 
880  min  of  initial  recovery  region.  These  control  regions  in  the 
time  domain  were  initialized  after  the  desired  Tc,max  =  42.4°C 
was  attained  in  the  hyperthermia  region  under  heat  stress.  The 
optimal  Ta  zone  was  followed  by  a  constant  Ta  =25°C  for  the 
rest  of  the  simulation  time.  A  GA  optimization  (10)  was  used 
to  minimize  a  given  cost  function  in  Eq.  11. 


Fig.  4.  Comparison  of  passive  cooling  rates  in  dead  mice.  The  plot  shown  here 
compares  the  experimental  cooling  rates  (solid  line)  with  the  computational 
mouse  model  (dashed  line).  Experimental  cooling  curves  of  4  different  mice 
are  plotted  around  the  mean  values.  Starting  with  the  same  initial  temperature 
of  42.7°C,  the  error  bar  represents  the  minimum  and  maximum  values  of 
temperature  in  different  mice  during  cooling. 


Fig.  5.  Distribution  of  temperature  across  a  mouse  body  under  heat  stress  at 
ambient  temperature  (Ta)  =  39.5°C.  Plot  corresponds  to  time  t  =  200  min 
under  hyperthermia  region.  A  minimum  temperature  of  15°C  is  observed  on 
the  front  legs  due  to  salivary  evaporative  heat  loss. 

T  T 

err+err+  ■  10  +  err_err _ 


where  V  is  the  value  of  the  cost  function,  err+  is  the  residual 
of  the  Tc  over  Tvet,  err-  is  the  residual  of  Tc  below  T.set,  and  N 
is  the  number  of  recorded  data  points  that  were  collected  from 
the  three-dimensional  solution.  The  range  of  Ta  that  can  be 
varied  during  the  optimization  is  restricted  to  (21°C  to  39°C). 
Five  GA  runs  were  conducted,  with  a  population  size  of  10, 
and  mutation  of  5%.  We  kept  members  of  previous  generations 
when  offspring  performance  was  inferior.  The  termination 
criteria  was  50  generations  without  change  in  the  value  of  the 
cost  function.  This  approach  was  presented  in  Grosman  and 
Lewin  (8).  The  optimized  parameters  were  represented  by  four 
binary  length  numbers,  which  lead  to  the  resolution  of  1.2°C. 

Figure  6A  depicts  a  comparison  of  thermoregulatory  profiles 
of  heat-stressed  mice  recovering  at  Ta  =  25°C  with  that  of 
different  Ta  values  obtained  by  GA  solution  and  other  subop- 
timal  cooling  profiles.  The  solution  obtained  by  the  GA  opti¬ 
mizer  (solid  line  A  and  B)  to  recover  the  Tt  back  to  Tjet 
suggested  that,  once  a  mouse  reaches  Tcanax,  it  should  recover 
at  Ta  =  27°C  for  110  min  initially.  In  the  second  stage,  the 
mouse  should  recover  in  a  mildly  heated  environment  of  Ta  = 
30°C  and  Ta  =  31°C  for  220  min,  followed  by  five  110-min 
segments  at  Ta  =  29,  28,  28,  27,  and  26°C,  representing  a 
gradual  reduction  of  Ta.  The  GA  solution  for  optimal  cooling 
reduced  the  Tc  after  heat  stress;  however,  once  the  Tc  dropped 
below  Tvet,  a  moderate  heating  profile  was  applied  to  prevent 
prolonged  hypothermia.  The  value  of  the  cost  function  pre¬ 
sented  by  Eq.  11  obtained  for  this  optimal  solution  was  3.8. 
The  dotted  lines  in  Fig.  6  (A  and  B)  describe  an  alternative 
solution  where  the  cooling  profile  started  with  initial  110  min 
at  Ta  =  21°C  and  then  Ta  =  25°C  for  the  rest  of  the  period. 
This  solution  resulted  in  a  longer  hypothermia  region  with  a 
cost  function  value  of  7.5.  The  dashed-dotted  line  in  Fig.  6  (A 
and  B)  describes  a  cooling  scenario,  where  the  mouse  was  kept 
at  Ta  =  25°C  throughout  the  recovery  period,  which  resulted  in 
slower  temperature  reduction  after  the  hyperthermia  and  a 
slower  recovery  from  the  hypothermia  with  the  cost  function 
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Fig.  6.  Optimal  Ta  values  obtained  by  genetic 
algorithm  solution  for  recovering  the  mice  Tc 
back  to  Tset  are  compared  with  suboptimal 
solutions  and  nominal  solution  after  the  max¬ 
imum  Tc  =  42.4°C  is  reached  in  240  min. 
Thermoregulatory  profiles  of  mice  exposed  to 
Ta  =  39.5°C  during  heat  stress  and  recovering 
at  Ta  =  25  °C  are  shown  in  A  with  dash-dot 
line,  and  the  corresponding  Ta  profile  is 
shown  in  B.  Mice  Tc  profile  using  optimal  Ta 
obtained  by  genetic  algorithm  solution  is 
shown  in  A  by  solid  line,  and  the  correspond¬ 
ing  Ta  controls  are  shown  in  B ,  which  are 
compared  with  suboptimal  solutions  in 
dashed  lines  in  both  plots. 


for  this  case  of  6.5.  In  the  last  scenario,  Ta  was  gradually 
reduced  from  Ta  =  33°C  to  25°C,  which  resulted  in  a  slower 
drop  from  hyperthermia  with  a  cost  function  of  10.2.  The  Tc 
recovery  time  for  a  heat-stressed  mouse  from  the  hypothermia 
region  with  constant  Ta  =  25 °C  was  around  12  h.  However, 
with  optimal  variable  Ta,  one  can  recover  the  mouse  Tc  back  to 
Tve!  in  around  7  h.  In  summary,  the  GA  solution  predicted  that, 
to  reduce  the  recovery  period,  the  mouse  should  be  cooled 
instantly  after  the  hyperthermia;  however,  after  the  mouse  Tc 
reaches  a  temperature  around  32°C,  a  moderate  heating  period 
that  gradually  is  reduced  to  Ta  =  25°C  should  be  applied. 

DISCUSSION 

Heat  exchange  with  the  environment  can  have  significant 
effects  on  thermal  tolerance  and  HS  recovery  in  mice.  A 
mammalian  thermoregulation  model  was  designed  to  simulate 
heat  stress  and  24-h  recovery  profiles  of  mice  with  consider¬ 
ation  of  passive  (conduction,  convection,  and  radiation)  and 
active  systems.  The  passive  system  was  developed  using  dif¬ 
ferent  thermophysical  and  physiological  properties  of  tissue 
materials,  and  the  active  system  based  on  a  set  point  integral 
controller  that  simulates  hypothalamic  control  of  thermoregu¬ 
lation.  Different  heat  transfer  mechanisms,  conduction,  con¬ 


vection,  radiation,  blood  perfusion,  evaporation,  and  metabolic 
heat  production,  are  considered  within  tissues  and  the  sur¬ 
roundings.  The  validity  of  the  passive  system  was  compared 
with  passive  cooling  profiles  of  mice  to  more  accurately 
simulate  mammalian  thermoregulatory  processes.  Using  Com- 
sol,  we  were  able  to  accurately  simulate  mouse  thermoregula¬ 
tory  profiles  during  heat  exposure  and  20  h  of  recovery.  The 
thermoregulation  model  of  the  mouse  shown  here  was  used  to 
reset  the  Tc  recovery  profiles  under  heat  stress  by  GA  opti¬ 
mizer. 

We  presented  a  novel  function  that  simulated  the  thermo¬ 
regulation  using  an  error  integrator  between  Tset  and  Tc.  The 
error  integration  action  that  manipulated  the  metabolic  activity 
was  motivated  both  by  early  observation  that  the  depth  and 
duration  of  the  hypothermia  is  directly  related  to  the  heat 
severity  (12),  and  to  maintain  an  offset-free  thermoregulation. 
Although  the  model  was  able  to  predict  the  depth  of  hypother¬ 
mia  and  recovery  times  fairly  well,  fever  during  the  second  day 
of  recovery  was  not  accurately  simulated.  Current  data  on 
post-heat  stress  Tc  elevations  suggest  that  fever  occurs  in 
response  to  endotoxin  leakage  across  damaged  gut  epithelial 
membranes  following  thermal  injury  (11).  Endotoxemia  causes 
inflammation  in  the  system  and  subsequent  production  of 
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cytokines,  which  are  known  to  mediate  thermoregulatory  re¬ 
sponses  (11).  The  morbidity  and  mortality  associated  with 
extreme  heat  stress  is  due  to  dysfunction  of  multiple  organ 
systems  as  a  consequence  of  the  systemic  inflammatory  re¬ 
sponse  syndrome  (1).  A  more  detailed  model,  including  these 
biochemical  changes,  should  be  able  to  predict  the  fever 
response  post-heat  stress. 

To  accurately  simulate  the  pre-heating  period,  Tset  had  to  be 
set  to  34°C,  and  then  changed  to  36.35°C  for  the  rest  of  the 
simulation.  The  change  in  T.ret  corresponds  to  the  thermoreg¬ 
ulatory  observation  in  small  rodents  and  is  successfully  pre¬ 
dicted  by  the  mathematical  model.  There  is  evidence  that  a 
hidden  set  point  of  a  unified  system  is  unnecessary,  and 
coordination  between  thermoeffectors  is  achieved  through  their 
common  controlled  variable,  body  temperature  (19).  However, 
in  the  present  work,  a  simplification  was  introduced  to  describe 
more  detailed  body  thermoregulation. 

Developing  effective  pharmacological  treatments  for  HS 
recovery  will  probably  require  relatively  demanding  multidis¬ 
ciplinary  research  of  the  biochemical  changes  occurring  at  the 
cellular  level.  On  the  other  hand,  examining  the  influence  of  Ta 
upon  the  recovery  of  animal  models  from  HS  conditions  is  a 
relatively  simpler  assignment.  Currently,  HS  is  being  treated 
by  rapid  cooling  therapy.  However,  while  a  rapid  reduction  of 
body  temperature  is  efficacious  following  HS  collapse,  when 
body  temperature  is  still  elevated,  efficacious  treatments  during 
ensuing  hours  and  days  of  recovery  have  yet  to  be  identified. 
The  relatively  detailed  three-dimensional  Comsol  simulation 
model  enabled  us  to  determine  the  efficacy  of  various  Ta 
treatment  strategies  for  HS  recovery.  A  GA  optimization  was 
conducted  to  minimize  the  offset  in  Tc  during  the  recovery 
period  using  eight  time  segments  of  110  min.  The  GA  optimi¬ 
zation  solution  predicts  that  using  Ta  =  27°C  for  the  initial  1 10 
min  of  recovery,  followed  by  Ta  =  30,  31,  29,  28,  28,  27,  and 
26°C  for  the  rest  of  the  segments,  resulted  in  a  7-h  hypothermia 
duration,  which  represented  a  reduction  of  58%  compared  with 
a  12-h  hypothermia  duration  at  a  constant  Ta  of  25°C.  The  in 
silico  model  developed  in  the  present  study  indentihes  a  series 
of  Ta  values  that  can  be  used  to  reestablish  T.set  following  HS. 
However,  hypothermia  may  be  a  thermoregulatory  response 
that  permits  survival  following  HS,  disruption  of  which  is 
detrimental  under  some  conditions  (12).  The  factor  10  in  Eq. 
11  was  chosen  to  emphasize  the  importance  of  reducing  the 
body  temperature  as  fast  as  possible.  There  is  no  obvious 
guidance  for  the  size  of  the  factor,  and  different  factors  could 
create  slightly  different  results.  However,  the  solution  that  was 
generated  suggests  rapid  cooling  after  the  heat  stress  condi¬ 
tions.  This  would  likely  be  generated  by  a  larger  factor  as  well. 

The  GA  is  a  stochastic  optimization  method  that  has  its  share 
of  drawbacks;  however,  understanding  the  nature  of  the  spe¬ 
cific  optimization  problem  and  the  lack  of  gradient-based 
information  lead  us  to  choose  the  GA  optimization  over  other 
alternatives.  A  gradient-based  optimization  would  lead  to  a 
large  number  of  evaluations  of  the  Comsol  model  that  would 
take  an  unreasonable  amount  of  time.  Moreover,  this  highly 
nonlinear  optimization  will  be  very  sensitive  to  the  initial 
conditions,  and  this  may  lead  to  a  yet  larger  number  of 
optimizations.  The  GA  provides  the  ability  to  conduct  a  search 
on  a  relative  rough  grid  of  temperatures  to  indicate  a  proof  of 
concept  showing  the  ability  of  a  simulation  to  predict  a  non¬ 
trivial  treatment. 


The  model  presented  is  a  simplification  of  reality.  The 
geometry  of  the  real  mouse  body  was  simplified.  The  meta¬ 
bolic  regulation  by  an  integrator  controller  mechanism  is  a 
hypothesis  that  could  incorporate  a  statistical  behavior  of  more 
complicated  regulation  mechanisms  in  reality.  We  based  the 
metabolic  heat  generation  on  Kleiber’s  law,  and  this  could  lead 
to  some  inaccuracies  when  used  across  species,  and,  indeed, 
we  needed  to  multiply  the  metabolic  heat  generation  by  a  factor 
of  three  to  enable  thermoregulation. 

All  mathematical  models  are  an  approximation  of  reality; 
therefore,  despite  the  many  advantages  of  using  these  models, 
one  has  to  take  into  account  the  limitations.  In  the  discussion, 
model  simplification  of  the  morphology  of  the  mouse  was 
done,  the  blood  system,  the  nervous  system  signals,  and  also 
the  heat  loss  through  the  respiratory  system  were  neglected. 
Moreover,  during  heat  stress,  small  rodents  adjust  body  tem¬ 
perature  on  a  hyperthermic  plateau.  The  height  of  the  plateau 
depends  on  salivation  threshold,  which  is  also  affected  by  Ta. 
This  mechanism  was  neglected  and  could  be  taken  into  con¬ 
sideration  in  future  versions  of  the  model.  Despite  numerous 
approximations,  the  model  is  enough  to  provide  accurate  pre¬ 
dictions  of  mouse  body  temperature  distribution  as  a  result  of 
changing  environmental  conditions. 

The  model  presented  here  demonstrates  the  usefulness  of  in 
silico  research  in  designing  novel  hypotheses  that  can  be  tested 
in  vivo.  Substantial  regarding  the  computer  model  in  terms  of 
HS,  in  silico  data  can  be  of  significant  benefit  in  narrowing  the 
identification  of  regulatory  windows  that  may  prove  efficacious 
in  treatment  of  HS  patients.  Testing  the  computer  model  via 
multiple  iterations,  complemented  by  the  inclusion  of  new  in 
vivo  data,  offers  a  powerful  experimental  tool  to  significantly 
enhance  the  discovery  of  novel  therapeutic  strategies.  Consis¬ 
tent  with  the  use  of  a  systems  biology  approach,  combining  in 
vitro,  in  vivo,  and  in  silico  models  holds  significant  promise  in 
reducing  animal  use,  as  well  as  the  ambiguity  associated  with 
experimental  protocols  consisting  of  a  large  number  of  vari¬ 
ables. 
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